Validation of IMU-based calculation


In [ ]:
import matplotlib
# matplotlib.use('Agg') #Comment this line out to see plots
import matplotlib.pyplot as plt
import numpy as np
import nvg.ximu.ximudata as ximudata
import nvg.ximu.compare_to_mocap as compare
import gc
#%matplotlib notebook #Uncomment to see plots

In [ ]:
from collections import defaultdict
from gc import get_objects
before = defaultdict(int)
after = defaultdict(int)
for i in get_objects():
    before[type(i)] += 1

In [ ]:
#import pdb
for ind in range(2,3):
    compare.main_2017(displacementComparison=True, orientationComparison=True,
                     trialindices=[ind],plotResults=False, anTime=60)

In [ ]:
for i in get_objects():
    after[type(i)] += 1

print [(k, after[k] - before[k]) for k in after if after[k] - before[k]]

In [ ]:
gc.collect()
plt.close("all")

In [ ]:
import sys
sys.getsizeof(compare.xdb)

In [ ]:


In [ ]:
import nvg.ximu.objgraph as objgraph
objgraph.show_most_common_types(limit=100)

In [ ]:
import matplotlib as mpl
mpl.style.available

In [ ]:
compare.xdb._pick_standing_reference("S12", "N")

In [ ]:
import nvg.io.qualisys_tsv as qtsv

In [ ]:
compare.xdb.estimate_joint_axes("S4", "N", ["LA", "LT"])

In [ ]:
compare.xdb.estimate_joint_center("S5", "N", ["LA", "LT"], startTime=220, anTime=40)

In [ ]:
import pdb
#compare.xdb.rotationEstimator = compare.xdb.GyroIntegratorOrientation(262.)
compdata = compare.markerdata_list()
compare.(compdata[-1:], imu="LA", markers=["ANKLE", "KNEE"], 
                  anTime=20, gThreshold=1e-1, lambda_gyro=1, 
                  lambda_incl=0.1,
                  var_angvel=1, var_incl=1,
                  resetAtIC=True, useCyclic=False, useLS=False)

In [ ]:
import pdb
#compare.xdb.rotationEstimator = compare.xdb.GyroIntegratorOrientation(262.)
compdata = compare.markerdata_list()
compare.main_2017(compdata[-1:], imu="LA", markers=["ANKLE", "KNEE"], 
                  anTime=20, gThreshold=1e-1, lambda_gyro=1, 
                  lambda_incl=1,
                  var_angvel=10, var_incl=1,
                  resetAtIC=False, useCyclic=True, useLS=False)

In [ ]:
pdb.pm()

In [ ]:
import pdb
#compare.xdb.rotationEstimator = compare.xdb.GyroIntegratorOrientation(262.)
compdata = compare.markerdata_list()
compare.main_2017(compdata[-1:], imu="LA", markers=["ANKLE", "KNEE"], 
                  anTime=20, gThreshold=1e-1, lambda_gyro=1, 
                  lambda_incl=0.1,
                  var_angvel=10, var_incl=1,
                  resetAtIC=True, useCyclic=True, useLS=True)

In [ ]:
import pdb
mdfiles = compare.markerdata_list()
compare.plot_marker_data(mdfiles[-1:], startTime=120, anTime=180)

In [ ]:
compare.xdb.hdfFile.close()

In [ ]:
#md = qtsv.loadQualisysTSVFile("/home/kjartan/Dropbox/projekt/nvg/data/solna09/S7/NVG_2012_S7_sync.tsv")

In [ ]:
#synctime = md.timeStamp + datetime.timedelta(seconds=5199.0/md.frameRate)

In [ ]:
a = (1,2,3,4)
a[-1:]
np.diff(a)

In [ ]:
import pdb
compare.main_2017()

In [ ]:
# Check the the correspondence between packet number and time
import numpy as np
from datetime import datetime, timedelta, date
fname = ("/media/ubuntu-15-10/home/kjartan/nvg/2012-09-19-S4/" 
        + "S4/LA-200/NVG_2012_S4_A_LA_00203_DateTime.csv")
packetTimeRaw = np.loadtxt(fname, dtype=np.int32, delimiter=',', skiprows=1)
times = [ datetime(year=pt_[1], month=pt_[2], day=pt_[3], 
                   hour=pt_[4], minute=pt_[5], second=pt_[6]) 
         for pt_ in packetTimeRaw ]
packetTime = np.asarray([ (packetTimeRaw[k,0]-packetTimeRaw[0,0], 
                           (times[k]-times[0]).total_seconds()) 
                         for k in range(len(times))])

In [ ]:


In [ ]:
Dp = (packetTime[-1,0] - packetTime[0,0])
Dt = (packetTime[-1,1] - packetTime[0,1])
freq = Dp/Dt
print freq

In [ ]:
from datetime import datetime, timedelta, date
timestampSync = datetime.strptime("2012-09-24, 16:52:59", "%Y-%m-%d, %H:%M:%S")
timestampRef = datetime.strptime("2012-09-24, 18:01:45", "%Y-%m-%d, %H:%M:%S")

floatSync = 267828.75835598
floatRef = 271955.08704553

dTime = timestampRef - timestampSync
dTimeFloat = floatRef - floatSync

print dTime
print dTimeFloat

In [ ]:
dTime.total_seconds()

In [ ]:
sfreq = compare.xdb.hdfFile.attrs['packetNumbersPerSecond']
print sfreq
type(sfreq)

In [ ]:
startTime = -120
anTime = 120
[la, t, s] = compare.xdb.get_imu_data("S12", "D", "LA", 
                                      startTime=startTime, anTime=anTime)
[ls, t, s] = compare.xdb.get_imu_data("S12", "D", "RA", 
                                      startTime=startTime, anTime=anTime)
accLA = la[:, 4:7]
accLT = ls[:,4:7]
plt.figure()
#plt.plot(np.mean(accLA**2, axis=1))
#plt.plot(np.mean(accLT**2, axis=1))
#plt.plot(accLA[:,-1])
#plt.plot(accLT[:,-1])
plt.subplot(211)
plt.plot(la[:,1:4])
plt.plot(ls[:,1:4])
plt.subplot(212)
plt.plot(la[:,4:7])
plt.plot(ls[:,4:7])

In [ ]:
gdir=np.mean(la[14000:140020, 4:7], axis=0)
print gdir
np.arccos(np.dot(gdir, [-1., 0, 0]))*180/np.pi

In [ ]:
syncLT = compare.xdb.get_PN_at_sync("S4","LT")
syncLA = compare.xdb.get_PN_at_sync("S4","LA")

print syncLT-syncLA
print ls[0,0] - la[0,0]

In [ ]:
import itertools
import pdb
[la, t, s] = compare.xdb.get_imu_data("S10", "N", "LA", split=True)
[lt, t, s] = compare.xdb.get_imu_data("S10", "N", "RA", split=True)
plt.figure()
for (la_, lt_) in itertools.izip(la[:40], lt[:40]):
    accLA = la_[:, 4:7]
    accLT = lt_[:,4:7]
    plt.subplot(211)
    #plt.plot(np.mean(accLA**2, axis=1))
    plt.plot(accLA[:,-1])
    plt.subplot(212)
    #plt.plot(np.mean(accLT**2, axis=1))
    plt.plot(accLT[:,-1])

In [ ]:
[la, t, s] = compare.xdb.get_imu_data("S4", "D", "LA", startTime=None, anTime=120)
la

In [ ]:
print la[0,0]
print syncLA

In [ ]:
[ra, t, s] = compare.xdb.get_imu_data("S4", "D", "RA")
accLA = la[:, 4:7]
accLT = ls[:,4:7]
plt.figure()
plt.plot(np.mean(accLA**2, axis=1))
plt.plot(np.mean(accLT**2, axis=1))

In [ ]:
ximudata.plot_sync()

In [ ]:
None < 0

In [ ]:
subj = compare.xdb.hdfFile['S4']
len(subj)

In [ ]:
imus = compare.xdb.hdfFile['S4']['N']

In [ ]:
imus

In [ ]:
for k in imus.iterkeys():
    print k

In [ ]:
a = np.arange(12)
a.shape=(4,3)

v = np.array([1,1,1])

np.cross(a, v)
print a

In [ ]:
print a[:-1,-1]
print a[:-1, -1][::-1]

In [ ]:
a = np.arange(12).reshape(4,3)
print np.linalg.norm(a, axis=1)
x = np.arange(3)
y = np.arange(3,6)
z = np.arange(6,9)
print np.stack( (x,y,z), axis=-1)

In [ ]:
np.reshape(np.insert(a, 0, [1,2,3]), (5, 3))

In [ ]:
a = np.arange(12)
ind = np.where(np.logical_and(a>4, a<9) )
a[ind]

In [ ]:
b = a.reshape(4,3)
np.dot(b, [1,1,1])

In [ ]:
b

In [ ]:
np.sum(b, axis=1)

In [ ]:
np.arctan2(Out[14], [0,0,0,0])

In [ ]:
b = {("Hej", "ha"):42, ("d","ha"):"s"}

In [ ]:
for (s_, t_) in b.keys():
    print s_
    print b[(s_, t_)]

In [ ]:
b.iteritems().next()

In [ ]:
for imu_ in compare.xdb.hdfFile["S4"]["N"]:
    print imu_

In [ ]:
s = "123,145.0, 150"
s1, s2, s3 = s.split(",")

In [ ]:
start = int(s2)
print start

In [ ]:
min(4, 14)

In [ ]:
a = np.arange(12)
a.shape=(4,3)
print a[:,0]
print a[0:4, 0]
print a[:,0:1]
a.shape=(12,)
print a[[0,-1]]

In [ ]:
from scipy.integrate import cumtrapz
b = a.astype(float)
print b
b.shape=(4,3)
print np.mean(b, axis=0)
b -= np.mean(b, axis=0)
print b
print cumtrapz(b, axis=0)

In [ ]:
a

In [ ]:
np.dot(a, np.arange(3))

In [ ]: